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Abstract 

We examine the basic mode structure of atomic Cooper pairs in an inhomogeneous Fermi gas. 
Based on the properties of Bogoliubov quasi-particle vacuum, the single particle density matrix 
and the anomalous density matrix share the same set of eigenfunctions. These eigenfunctions 
correspond to natural pairing orbits associated with the BCS ground state. We investigate these 
orbits for a Fermi gas in a spherical harmonic trap, and construct the wave function of a Cooper 
pair in the form of Schmidt decomposition. The issue of spatial quantum entanglement between 
constituent atoms in a pair is addressed. 

PACS numbers: 03.75.Ss, 05.30.Fk, 74.20.Fg, 03.67.Mn 
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I. INTRODUCTION 



Quantum degenerate gas of two-component Fermi atoms provides a well-controlled sys- 
;em for studying Cooper pairing responsible for superfluidity phenomena in the BCS regime 
With the advances in cooling and trapping_^ techniques, recent experiments have 
demonstrated various consequences of pairing ji! 0, 111 0] . O ne of the major tasks for theo- 
retical investigation is the structure of Cooper pairs in non-uniform finite systems. Since the 
weakly interacting ultra-cold atoms mainly involve two-body s-wave scattering processes, de- 
tailed analysis can be performed via the pseudo-potential method and Bogoliubov-deGennes 
(BdG) equations jz,!^- An important quantity that can be obtained self-consistently by this 
approach is the gap function. Such a quantity measures the pairing field at a given point 
in space, but it lacks the information about correlation between atoms at different spatial 
points. In order to gain a more complete picture of pairing, it is useful to examine properties 
of two-point correlation functions 9|. 

In this paper we examine the underlying mode structures inherent in two-point correlation 
functions: (r^^^)) and (^ / g(r 1 ) , Q ,(r2)), with ipj (j = being the field operator 

associated with the spin component j. For symmetric systems, we employ the fact that these 
correlation functions are built up by the same set of orthogonal eigenfunctions [ljj]. Such 
eigenfunctions can be interpreted as natural pairing orbits that form the BCS wave functions 
in inhomogeneous systems. We note that the standard textbook description of BCS wave 
function refers to infinite homogeneous systems such that each pairing orbit corresponds to 
the eigenfunctions of opposite momenta [8] . However, the presence of a confining potential in 
inhomogeneous systems could drastically alter the pairing orbits. Although the time-reversal 
symmetry helps to fix a certain set of quantum numbers in pairs, the exact form of pairing 
orbits are difficult to find in general. One possible solution is to treat the pairing orbits 
as (unknown) variational functions. With the usual BCS ansatz and variation technique, a 
set of nonlinearly coupled equations of pairing modes have been derived Q, 11 ] . However, 
these nonlinear equations have no simple physical in terp retations, and numerical solutions 
have only been demonstrated in nuclear systems |l(3, As the particle numbers in Fermi 
gases systems are typically much larger than that in nuclei, numerical method based on 
variational method could become difficult. To our knowledge, natural orbits of Cooper pairs 
in a trapped Fermi atomic gases have not been fully explored. 
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One of the main purposes of our paper is to indicate an efficient way of determining 
pairing orbits, directly from the two-point correlation functions obtained by BdG mean field 
equations. To illustrate our method, we will examine the gas in a spherical harmonic trap. 
For such systems, Bruun and Heiselberg have provided useful insight about the approximate 
forms of pairing ««, with a d^er/approach QQ. Here, « wiH present —1, 
exact examples of pair orbits. With these pairing orbits we can further construct and study 
the spatial wave function of a Cooper pair. In addition, we will address the quantum 
entanglement between two constituent atoms in a pair. Recently, one of us have addressed 
the issue of quantum entanglement of two atoms due to s-wave scattering [1,41]. The study 
of a Cooper pair would shed some light on the importance of many-body effects. 

II. THEORY 

A. The model and BdG equations 

To begin with, we write down the model Hamiltonian describing an interacting two- 
component Fermi gas in a trapping potential Uo(r), 



where Hq = — <^V 2 + Uq(t) — fi, with m and fi being the particle mass and chemical 
potential respectively. The coupling strength g is related to the s-wave scattering length a 
via g = AiTh 2 a/m. In this paper we assume that a is negative and the number of atoms in 
each component is the same. 

Under the mean field approximation, the Hamiltonian can be diagonalized through the 
Bogoliubov transformation: ip a (r) = Y, v u n{ r )lm ~ ^OOtJ/jj and ^( r ) = Y, n u v( r hvP + 
u*(r)7^ a , such that the ground state is the vacuum state of the quasi-particle operators 
7,,'s, i.e., Bogoliubov vacuum. The quasi-particle wave functions u v (r) and v v (r) are solved 
self-consistently by the Bogoliubov-deGennes (BdG) equations j^: 



H 




(1) 



[F (r) + W(r)} u v {r) + A(r)u„(r) = E n u v (r) 
A*(r)u v ( r )-[H (r) + W(r)]v v (r) = E v v v (r) 



(2) 
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with W(r) = g'Yl.n l^fr)! 2 , and A(r) = — g J2 V u v( r ) v ri( r )- ^ ^ s important to note that A(r) 
has a 1/r divergence m Ujl. Such a divergence can be eliminated by replacing q with a 

nn 

regularized effective coupling constant as described in Ref . |15|, |lg] . 



B. Two-point correlation functions 

Given that the ground state of the system is described by the Bogoliubov quasi-particle 
vacuum, the normal density matrix p(ri,r 2 ) = (^(ri)^^)) = (^Uri)^^)) and the 
anomalous density matrix u{yi,y 2 ) = (ipp(ri)ip a ( y r 2 )) take the form: 

p(r u r 2 ) = ^ v n( T i) v n( T 2) ( 3 ) 
v(r u r 2 ) = - ^ ^(ri)^(r 2 ). (4) 

In this paper we assume that p(r l5 r 2 ) and z/(r!,r 2 ) are real symmetric matrices. Such a 
symmetric property can be shown explicitly in spherically trapped systems that we will 
discussed later in the paper. 

From the properties of quasi-particle wave functions, it can be shown that p and v com- 
mute (see Appendix), i.e., 

J d 3 r 1 [p(r,ri)z/(ri,r 2 ) - z/(r, ri)p(ri,r 2 )] = 0. (5) 

This implies that the normal density matrix and the anomalous density matrix share a 
common set of eigenfunctions defined by the integral equations: 

J d 3 r 2 p(ri,r 2 ) / n (r 2 ) = Kfn^i) (6) 

J d 3 r 2 i/(n, r 2 ) / n (r 2 ) = Xnfn{n) (7) 

where A n and Xn are real eigenvalues, and {/„} is a complete and orthogonal set of eigen- 
functions. Furthermore, p(ri,r 2 ) and z/(ri,r 2 ) can be expressed as the following bilinear 
expansion: 

(V4(ri)Var 2 )) = E n A «/n( r i)/:(r 2 ) (8) 
(M^M^)) = Xnfn(ri)f*(r 2 ). (9) 

Here the convergence is ensured by the squared integrable property of the kernels, according 
to the theory of integral equations. 
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C. Natural orbits 



The eigenfunctions f n are called natural orbits. To interpret such orbits, let us write 
down a general mode expansion of operators ip a and ipp such that: ip a {r) = ^„ (r)a n , 
^(r) = <fin\r){3 n , where {0^} and {4>n } are two sets of complete orthogonal functions 
to be determined, and a n and j3 n are the corresponding fermion annihilation operators. Next 
we construct the BCS ground state in the standard form: 

\<$>) = \[{u n + v n aiPi)\V) (10) 

n 

with the normalized coefficients u n and v n satisfying \u n \ + \v n \ =1. From the fact that 

Wt(ri)V>«(r 8 )|<&) = ^(ri)0i Q) (r 2 ) \v n \ 2 (11) 
($|^(r 1 )^(r 2 )|$) = ^ € ] {vxW{v 2 )u* n v n (12) 

we immediately see that the constructed BCS ground state (fTU|) is consistent with the BdG 
results Eq. © and © if 

<Pt\r) = f:(r) (13) 
0i /3) W = /n(r) (14) 

and v n = \/\ n and u n = Xn/V^n (Appendix). In other words, /„ and its conjugate are 
indeed the pairing modes needed for the construction of the BCS ground state. This is 
shown by matching the correlation functions obtained from BdG mean-field equations. In 
fact, with the help of Wick's theorem, the choice of mode functions (13) and (14) also match 
higher order correlation functions. Hence the equivalence between Bogoliubov vacuum and 
the BCS ground state in a trapped system can be established explicitly. 

Historically, the concept of natural orbits were applied to pairing problems of nucleons 
and one can find some early references in Ref . |ipj . Here we employ such an idea to analyze 
the structure of Cooper pairs in atomic systems. We should emphasize that the key to 
the existence of /„ is the symmetric property of p, which makes p and v commute. In 
the appendix we indicate this point in a simple derivation of Eq. (j3J). For asymmetric 
systems, such as the systems with imbalance populations of the two components, p and v 
do not commute in general. This would then forbid the ground state to be in the BCS form, 
although the BdG equations could still be used to describe asymmetric systems. 
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D. Cooper pair wave function and quantum entanglement 



The description of a Cooper pair is inherited from the mean field description, where every 
pair in the system is assumed to be identical. Let us now construct the wave function of 
a Cooper pair based on the pairing orbits. This is achieved by noting that the BCS state 
takes the form: 

|$) oc exp(5>,a}/3j)|0> = EtT |0) (15) 

j k=0 

where Kj = Vj/v,j, = J2j K j a ]Pj, an d j stands for quantum numbers collectively. The 
wave function corresponding to the pair creation operator is given by 

F(r 1 ,r 2 ) = C^^/ J (r 1 )/;(r 2 ) (16) 

j 

where C is a normalization constant. Such a function is interpreted as the wave function of a 
Cooper pair. One reasoning is based on the fact that the BCS state is dominated by particle 
numbers near the mean value N when the particle number is large, i.e., k w N terms in Eq. 
(15) contribute most. This allows us to have an approximate picture of N pairs described 
by A^lO). It is worth noting that each term in expansion (16) is weighted by the factor Kj, 
which is different from Eq. (8) and (9). Since Kj is larger for the more filled orbits, orbits 
below the Fermi level could have more contributions individually. 

Eq. (16) is precisely a form of Schmidt decomposition of a bipartite system This 
is observed by the fact that fj are orthogonal basis functions, and Kj are corresponding 
Schmidt eigenvalues. Schmidt decomposition provides a useful way to characterize quantum 
entanglement between two subsystems in pure states. In particular the degree of quantum 
entanglement can be learned by the von Neumann entropy according to the distribution of 
Schmidt eigenvalues. Here we apply this method to address the entanglement between two 
atoms in a Cooper pair. A transparent and direct measure of entanglement is the 'average' 
number of Schmidt modes involved. The effective Schmidt number K provides this average 



14, 



19J: 



^ = 1/£K! 4 , (17) 

3 

where Kj = Ckj. Such a quantity is also the inverse of the purity function. The larger the 
value of K, the higher the entanglement. Note that the entanglement here refers to the 
spatial degrees of freedom, which should be distinguished from the spin entanglement [l8|. 



III. FERMI GASES IN A SPHERICAL HARMONIC POTENTIAL 



Having described the formalism, we now proceed to study natural orbits of a Fermi gas 
confined in a spherical harmonic potential: Uq(t) = |ma; 2 r 2 , with uj being the trapping fre- 
quency. First, we determine the quasi-particle mode functions u v 7 s and v^s numerically from 
the BdG equations. The spherical symmetry of the system gives: ru v i m (r) = u v i(r)Yi m (9 , <f>), 
r ?Vm( r ) = v n i{r)Yi m {9 , 0), with Yi m being the spherical harmonic functions. It should be 
noted that in order to remove the ultraviolet divergence, we employ an effective coupling 
constant g c s{v), such that the whole set of equations is independent on cutoff. The formu- 
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16]. 



and we will skip 



lation of renormalization scheme has been discussed extensively in 

the details in this paper. 

Once numerical solutions of u^s and i^'s are found, the correlation functions p(r 1; r 2 ) 

and z/(i7,r 2 ) are obtained from Eqs. Q and A significant simplification can be made 

i 

by using the addition formula: {21 + l)P/(cos7) = 4n ^ (#i, 0i) Y/ m (6* 2 , 2 ), where 

m=— I 

Pi (cos 7) is the Legendre polynomial, and 7 is the angle between between ri and r 2 . This 
gives 

(^(rOVafa)) = ^P(cos 7 )a*(r-i,r 2 ) (18) 
<^( ri )^(r 2 )> = -^P(cos 7 )A(ri,r 2 ) (19) 

where 

aj(ri,r 2 ) = — > '- 20 

Att r x r 2 

2/ + l ^ ^(n)^(r 2 ) 

A(ri,r 2 ) = — — > • 21 

Air ' rir 2 
v 

Therefore the correlation functions Q and (J3J) depends only on three variables: 7*1, r 2 and 7. 
The integral eigenvalue equations © and (J7J), can be solved conveniently by expressing the 
eigenfunctions as / n zm( r ) = ^h nl (r)Yi m (9, 0), with h n i(r) being the radial function associated 
with the radial quantum number n and orbital angular momentum I. Equations (6) and (7) 
then become, 

r\r 2 ot l {r u r 2 )h n i{r 2 )dr 2 = |^| 2 /i rai ( r i) ( 22 ) 
rir 2 (3i(ri,r 2 )h nl (r 2 )dr 2 = u* nl v n h n i(ri) (23) 




FIG. 1: (Color online) The g= -1.15 (solid line) and the non-interacting (dashed line) radial 
density distribution function, normalized to N= f 4irr 2 p(r)dr. Inset shows the pairing field A(r) 
with g= -1.15. 

which are one-dimensional integral equations, as only the radial coordinates are involved. 

To provide a concrete example, we consider a system of N = N a + Np = 2N a = 10912 
particles with the interaction strength g = —1.15 in trap units, (i.e., energy in fuv, length in 
y^h/muj). The choice of such a particle number corresponds to the Fermi energy ep = 31.5Hui 
in the non-interacting limit. Eq.Q was solved in a truncated harmonic oscillator states 
basis, with a cutoff energy ~ 180/iu; £L_To achieve convergent results, we employed the 
regularization scheme according to Ref. The particle density p(r) and pairing potential 
A(r)are consistently solved and shown in FigGJ In this example, A(0) = 7.2huj represents 
a modest strong coupling. Comparing with the non- interacting system (dashed line), we 
see that the particles are significantly dragged towards the center of the trap due to the 
attractive interaction. 

After solving the integral eigenvalue equations for h nh we obtain the distributions of 

eigenvalues \v n i\ 2 and u* nl v n i for various angular momentum quantum number I. We note 

that in the case of non-interacting systems at zero temperature, is a step function, i.e., 

|^nz| 2 = ©(A 4 — | ~ I — 2n). This corresponds to the fact that non-interacting atoms fill up 

all the trap energy states up to the Fermi level. For the interacting system considered here, 

the sharp edge of the step function is smeared out as shown in Fig. 2a. Such a smearing 

(1 

effect is similar to what appears in uniform BCS systems [8J. However, the difference here 
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FIG. 2: (Color online) Occupation number \vki\ 2 (upper panel) and the pairing amplitude u n iv n i 
(lower panel) with 1=0 (diamond), 5 (triangle), 20 (circle). Same parameters as in Fig. 1. 



is that /„/ are now the pairing basis instead of plane waves. 

In Fig. 2b, we show the quantity u^fini, which measures the coherence between occupied 
and un-occupied states. We see that u* nl v n i reaches a peak at a certain values of radial 
quantum number n. These n's correspond to orbits f n i that have average energies close to 
the chemical potential, and hence for higher /, the peak appears at lower n. The orbits near 
the peak contribute most significantly to the gap function. The shapes of some of these 
orbits are plotted in Fig. El where the radial part \h n i(r)\ 2 at various quantum numbers 
I is shown. Comparing with bare eigenfunctions of the trap (dashed line) with the same 
quantum numbers n and /, we observe similar oscillatory patterns but the envelopes are more 
concentrated towards the trap center. This indicates that the trap's eigenfunctions do not 
provide a good approximation to the actual pairing orbits, at least for the moderate strong 
coupling considered here. However, for a much weaker coupling A(0) <C hw (not shown), 
we do find a good agreement between pairing orbits and the bare trap's eigenfunctions near 
the Fermi surface, which is expected according to the argument in Ref . 12 1 . 

In Fig. HJ we illustrate the shape of the Cooper pair wave function by plotting the 
quantity P = |rir2-F(ri, r2, 7)| 2 at various angular separation 7. Note that P is proportional 
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FIG. 3: (Color online) Some radial mode functions \h n i(r)\ 2 of natural orbits, normalized to 
/ \h n i{ r )\ 2 dr = 1, with g= -1.15 (solid line). These functions correspond to those pairing orbits 
with a peak value of u n \v n i ~ 0.5 at a given I. The dashed lines correspond to the bare harmonic 
trap radial wave functions having the respective quantum numbers. 

to probability density of simultaneously finding the two particles at r x and r 2 , and the 
weighting factor rir 2 is included for spherical systems. In both Fig. 4a and Fig. 4b, we see 
some interesting fringes patterns, but the main feature is that the function P appears to 
be localized near the diagonal r\ = r 2 , indicating that both atoms are likely to be found in 
the same radial distance from the trap center. In addition, by comparing Fig. 4a and 4b, 
P drops significantly when angular separation 7 increases. This suggests that F(ri,r 2 ,'j) 
mainly concentrates at = r 2 , not just the same radial distance. For the example shown 
in Fig. 4a, the width of P near the peak is about 1.2 which is comparable to the bare trap 
ground state. 

The effective Schmidt number K defined in Eq. (17) is found to be K « 1749, which 
is roughly the same order of particle number in the trap of this example. It is useful to 
compare this number with the value of K of a two-atom system with the same scattering 



length. According to the calculation in Ref. 1J], a two-atom system in the ground state has 



a very weak entanglement if the scattering length is small compared with the trap length 
unit. For the parameters used in Fig. 4, the scattering length is a = —0.09 which gives 
K « 1 in the two-atom system. Therefore a Cooper pair, which is due to many-body effects, 
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FIG. 4: (Color online) An illustration of the spatial profile of a Cooper pair, where 
P = \ r i r 2F( r i, r 2,7)\ 2 1 with (a) 7 = and (b) 7 = tt/20. Same parameters as in Fig. 1. 



processes a much stronger entanglement. 



As a final remark, we note in Ref. 20] that the ratio K/Nj (where N a = Np is particle 
number of a component) is suggested to be an indicator of determining whether a composite 
two-particle system behave as a boson or not. Specifically, K/N ^> 1 is the regime where 
composite particles can be described by creation and annihilation operators obeying bosonic 
commutation relations. Here in our example, we have K/N pa 1749/5456 which is still 
smaller than one. Therefore the entanglement is not strong enough to hold the Cooper pair 
together boson. 



IV. CONCLUSIONS 



To summarize, eigen-mode expansion of correlation functions provides a powerful tool to 
reveal the underlying coherent structures. Such a strategy is also known useful in various 
areas of physics, such as nuclear physics and optics [21]. In this paper, we have examined 
the eigenf unctions of two-point correlation functions (3) and (4) that can be obtained from 
the solutions of BdG equations. These eigenfunctions correspond to natural pairing orbits in 
the BCS state (10), and they serve as Schmidt basis vectors for the construction of Cooper 
pair wave functions. Hence, one can further analyze the quantum entanglement associated 
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with the spatial degree of freedom. We have demonstrated the method in a spherically 
trapped system, in which the natural orbits are calculated explicitly. In particular, our 
numerical results indicate features that reflect the strong quantum entanglement between 
two constituent atoms in a pair. It is useful to point out that the method could also 
be extended to BEC-BCS crossover problems, where BdG equations are modified by the 
molecular field and a hybrid form of BCS-molecule state may be constructed. The structure 
of natural orbits in this regime is a very interesting topic for future investigations. 
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APPENDIX 

The vanishing commutator given in Eq. (jSJ) can be derived under a rather general con- 
sideration (see for example, in Ref. Q]). Here we re-derive the result in a more transparent 
way. We start with the orthogonality and completeness relations of the quasi-particle wave 
functions M 




(24) 



(25) 



Together with Eqs. (j3J) and (JU), we can show 




mn 




run 




mn 




(26) 
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Now we see that Eq. (jSJ) holds if the symmetric condition: p(ri,r2) = p(r2,ri) is satisfied. 
With Eqs. (}2"4"j) and (J23j) and the symmetric p(ri,r2), we also have the relation: 



Projecting both sides of Eq. (|27jl onto /n(r 2 ), we obtain a relation of the eigenvalues, 



which is consistent with the normalization condition \u n \ 2 + |t>J 2 = 1 of BCS wave function. 
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